Quasi-adiabatic quantum Monte Carlo algorithm for quantum evolution in imaginary time 



o 

(N 
O 

Q 



Cheng-Wei Liu, Anatoli Polkovnikov, and Anders W. Sandvik 

Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA 

(Dated: December 20, 2012) 

We propose a quantum Monte Carlo (QMC) algorithm for non-equilibrium dynamics in a system with a 
parameter varying as a function of time. The method is based on successive applications of an evolving Hamil- 
tonian to an initial state and delivers results for a whole range of the tuning parameter in a single run, allowing for 
access to both static and dynamic properties of the system. This approach reduces to the standard Schrodinger 
dynamics in imaginary time for quasi-adiabatic evolutions, i.e., including the leading non-adiabatic correction 
to the adiabatic limit. We here demonstrate this quasi-adiabatic QMC (QAQMC) method for linear ramps of the 
transverse-field Ising model across its quantum-critical point in one and two dimensions. The critical behavior 
can be described by generalized dynamic scaling. For the two-dimensional square-lattice system we use the 
method to obtain a high-precision estimate of the quantum-critical point (h/J) c = 3.04463(12), where h is the 
transverse magnetic field and J the nearest-neighbor Ising coupling. The QAQMC method can also be used to 
extract the Berry curvature and the metric tensor. 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods 1 - 2 have become in- 
dispensable tools for ground-state and finite-temperature stud- 
ies of many classes of interacting quantum systems, in par- 
ticular those for which the infamous "sign problem" can be 
circumvented. In ground-state projector methods, an op- 
erator P(f5) is applied to a "trial state" \^o), sucn that 
\^>p) = P(y8)|^o) approaches the ground state of the Hamil- 
tonian H when f3 — > oo and an expectation value (A) = 

p\A\$> p) / Z , with the norm Z = (j&p\iflp), approaches 
its true ground-state value, (A) — > (0|A|0). For the pro- 
jector, one can use P(j3) = exp (—f3H) or a high power of 
the Hamiltonian 3 , P(M) = (— TL) M . Here we will discuss 
a modification of the latter projector for studies of dynamical 
properties of systems out of equilibrium. 

Our work is a further development of the method introduced 
recently in Ref. [4], where it was realized that a modification 
of the projector approach with P(f3) = cxp (~f3H) can be 
used to study non-equilibrium set-ups in quantum quenches, 
where a parameter of the Hamiltonian depends on time ac- 
cording to an arbitrary protocol. By performing a standard 
Wick rotation of the time axis, a wave function is governed by 
the Shrodinger equation in imaginary time t = —ir (r being 
real), 

d T \4ir)) = -H[X(t)Mt)). (1) 

Here the Hamiltonian depends on the parameter A through 
time, e.g., 



H = H + A(t)V, 



(2) 



where V and Hq typically do not commute. The method is not 
limited to this form, however, and any evolution of % can be 
considered. The Schrodinger equation has the formal solution 

|V(r)) = C/(t)|V>(t )), (3) 

where the imaginary-time evolution operator is given by 



where T T indicates time ordering. A time-evolved state 
U(t)\^(tq)) and associated expectation values can be sam- 
pled using a generalized projector QMC algorithm. In 
Ref. [4] it was demonstrated that this non-equilibrium QMC 
(NEQMC) approach can be applied to study dynamic scal- 
ing at quantum phase transitions, and there are many other 
potential applications as well, e.g., when going beyond stud- 
ies of finite-size gaps in "glassy" quantum dynamics and the 
quantum-adiabatic paradigm for quantum computing. 

Here we introduce a different approach to QMC studies of 
quantum quenches which gives results for a whole range of 
parameters A G [X(tq), A(t)] in a single run (instead of just 
the final time), at a computational effort comparable to the 
previous approach. Instead of using the conventional time- 
evolution operator Eq. (4), we consider a generalization of the 
equilibrium QMC scheme based on projection with (— H) M , 
acting on the initial ground state of %[A(to)] with a product 
of evolving hamiltonians; 

Pali = [-«(A M )]....[-«(A a )][-«(Ai)], (5) 



where 



A( = Aq + tA\, (6) 



and A a = [At+i — Xt)/M is the single-step change in the 
tuning parameter. 5 Here we will consider a case where the 
ground state |$(Ao)) of H(Xq) is known and easy to generate 
(stochastically or otherwise) and the ground state for other A- 
values of interest are non-trivial. The stochastic sampling used 
to compute the evolution then takes place in a space represent- 
ing path-integral-like terms contributing to the matrix element 
(the norm) (*(Ao)|Pi,m-Pm,i|*(Ao)). We will also later con- 
sider a modification of the method in which the ground state 
at the final point A m is known as well, in which case contri- 
butions to (*& (Xm)\Pm, i\^(Xo)} are sampled. 

Staying with the doubly-evolved situation for now, we eval- 
uate generalized expectation values after t out of the M oper- 
ators in the product (5) have acted; 



U(t) = IV exp 



dr'nlXir')} 



(4) 



(A)t 



OP(Ao)|Pi,A/P A /, f+ iAP t , 1 |*(A )) 

(*(A )|Pi,MiVi|*(A )) 



(7) 
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We will refer to this matrix element as an asymmetric expec- 
tation value, with the special case t = M corresponding to a 
true quantum mechanical expectation value taken with respect 
to an evolved wave function 



Pm,i|*(A )) 



v / (*(Ao)|A,mPm,i|*(Ao))' 



(8) 



which approaches the ground state |^[A(tjv/)]) of the Hamil- 
tonian H[X(tm)} for M — > oo. 

As we demonstrate in detail in Sec. II, a quench veloc- 
ity v oc A^A can be defined such that the symmetric ex- 
pectation valye (A) t= M in Eq. (7) approaches the expecta- 
tion value (A(t = t)) after a conventional linear imaginary- 
time quantum quench with Eq. (4) done with the same ve- 
locity v, if v is low enough. In fact, the two quantities are 
the same to leading (linear) order in v, not only in the strict 
adiabatic limit v — > 0. We therefore name this scheme the 
quasi-adiabatic QMC (QAQMC) algorithm. Importantly, the 
leading corrections to the adiabatic evolution of the asymmet- 
ric expectation values for any t contain important information 
about non-equal time correlation functions, very similar to the 
imaginary-time evolution. 

The principal advantage of QAQMC over the NEQMC ap- 
proach is that expectation values of diagonal operators in the 
basis used can be obtained simultaneously for the whole evo- 
lution path Ao . . • Am, by measuring (A) t in Eq. (7) at arbi- 
trary t points. The QAQMC scheme is also easier to imple- 
ment in practice than the NEQMC method because there are 
no time integrals to sample. 

As mentioned above, we will here have in mind a situation 
where the initial state |\l/(Ao)) is in some practical sense "sim- 
ple", but this is not necessary for the method to work — any 
state that can be simulated with standard equilibrium QMC 
methods can be used as the initial state for the dynamical evo- 
lution. The final evolved state \ipAi) can be very complex, 
e.g., for a system in the vicinity of a quantum-critical point or 
in a "quantum glass" (loosely speaking, a system with slow 
intrinsic dynamics due to spatial disorder and frustration ef- 
fects). Here, as a demonstration of the correctness and util- 
ity of the QAQMC approach, we study generalized dynamic 
scaling in the neighborhood of the quantum phase transitions 
in the standard one-dimensional (ID) and 2D transverse-field 
Ising models (TFIMs). 

As noted first in Ref. [4] the NEQMC method can be used 
to extract the components of the quantum metric tensor, 6 the 
diagonal elements of which are the more familiar fidelity sus- 
ceptibilities. Thanks to its ability to capture the leading non- 
adiabatic corrections to physical observables, the QAQMC 
approach can also be used for this purpose, and, as we will 
discuss briefly here and in more detail in Ref. [7], one can 
also extract the Berry curvature through the imaginary anti- 
symmetric components of the geometric tensor 

The rest of the paper is organized in the following way: In 
Sec. II we use adiabatic perturbation theory (APT) to demon- 
strate the ability of the QAQMC scheme to correctly cap- 
ture the standard Schrodinger evolution in imaginary time, not 
only in the adiabatic limit but also including the leading cor- 
rections in the quench velocity. We show how these leading 



corrections correspond to the geometric tensor. In Sec. Ill we 
discuss tests of the QAQMC scheme on ID and 2D TFIMs, 
and also present a high-precision result for the critical field in 
the 2D model. In Sec. IV we summarize our main conclusions 
and discuss future potential applications of the algorithm. 



II. ADIABATIC PERTURBATION THEORY 

The key question we address in this section is whether the 
matrix element (A) t in Eq. (7) can give useful dynamical in- 
formation for arbitrary "time" points t in the sequence of 2M 
operators. The expression only reduces to a conventional ex- 
pectation value at the symmetric point t = M, and even there 
it is not clear from the outset how (A) t= M computed for dif- 
ferent M relates to the velocity dependence of the expecta- 
tion value {^(0)\U*(t)U(t)\^(0)) based on the Schrodinger 
time-evolution operator in Eq. (4). Going away from the sym- 
metric point brings in further issues to be addressed. For in- 
stance, there is no variational property of the asymmetric ex- 
pectation value (H)t of the Hamiltonian for t ^ M. Nev- 
ertheless the approach to the adiabatic limit is well behaved 
and we can associate the leading deviations from adiabatic- 
ity with well defined dynamical correlation functions which 
appear as physical response in real time protocols. We show 
here, for the linear evolution Eq. (6), that one can identify a 
velocity v oc N/M such that a linear imaginary-time quench 
with A t = vt in Eq. (6) gives the same results in the two ap- 
proaches when t = M, including the leading (linear) correc- 
tions in v. For t ^ M, the relevant susceptibilities in QAQMC 
defining non-adiabatic response are different than at t = M 
but still well defined and obeying generic scaling properties. 

In order to facilitate the discussion of the QAQMC, we here 
first review the previous APT approach for Schrodinger imag- 
inary time dynamics 4 - 7 and then derive analogous expressions 
for the product-evolution. After this we discuss some proper- 
ties of the symmetric and asymmetric expectation values. 



A. Non-equilibrium Quantum Monte Carlo 

The NEQMC method 4 uses a path-integral-like Monte 
Carlo sampling to solve the imaginary-time Shcrodinger equa- 
tion Eq. (1) for a Hamiltonian H[A(t)] with a time depen- 
dent coupling. The formal solution at time r is given by the 
evolution operator Eq. (4). In the strict adiabatic limit the 
system will follow the instantaneous ground state, while in 
the slow limit one can anticipate deviations from adiabatic- 
ity which will become more severe in gapless systems and 
in particular near phase transitions. Let us discuss the lead- 
ing non-adiabatic correction to this imaginary-time evolution. 
The natural way to address this question is to use APT, similar 
to that developed in Refs. [8] and [9] in real time. We here fol- 
low closely the discussion of the generalization to imaginary 
time in Ref. [4]. 

We first write the wave function in the instantaneous eigen- 
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basis {|n(A))} of the time-dependent Hamiltonian H[X(t)]: 
^(t)=]>>„(t)KA(t))>. (9) 

We then substitute this expansion into Eq. (1), 

da ri 



dr 



^o m (r)(n|9 T |m) = -f„(A)a„(r), (10) 



where £ n (X) are the eigenenergies of the Hamiltonian H(\) 
corresponding to the states |n) for this value of A. Making the 
transformation 



a„(r) = a„(r) exp J £ n (r')dr' 
we can rewrite Eq. (1) as an integral equation; 

a n (r) = a n (0) + / dr' {n\d T ,\m)a m (T') 



(ID 



x exp 



dr" (£ n (r") - £ m ( T ") 



, (12) 



where it should be noted that a„(0) = a n (0). In principle we 
should supply this equation with initial conditions at r = To, 
but this is not necessary if |t | is sufficiently large, since the 
sensitivity to the initial condition will then be exponentially 
suppressed. Instead we can impose the asymptotic condition 
a„(r — s- — oo) — > S n o, which implies that in the distant past 
the system was in its ground state. 

Eq. (12) is ideally suited for an analysis with the APT. In 
particular, if the rate of change is very small, A(r) — > 0, then 
to leading order in A the system remains in its ground state; 
a m (r) ps S m o (except during the initial transient, which is not 
important because we are interested in large | To [ ) . In the next 
higher order the transition amplitudes to the states n ^ are 
given by; 



a n (0) ps — / dr (n|9 T |0) exp 



- / dT'A n0 (T') 



(13) 



where A„o(t) = £ n {r) — £o(t). The matrix element above 
for non-degenerate states can also be written as 



(n\d T \0) = -(n\d T U(T)\0}/A n0 {T). 



(14) 



In what follows we will assume that we are dealing with a 
non-degenerate ground state. 

To make further progress in analyzing the transition am- 
plitudes Eq. (13), we consider the very slow asymptotic limit 
A — >• 0. To be specific, we assume that near t = the tuning 
parameter has the form (see also Ref. [9]); 



A(t) ps A(0) + 



v\\t t \ 



(15) 



The parameter v\, which controls the adiabaticity, plays the 
role of the quench amplitude if r = 0, the velocity for r = 1, 



the acceleration for r = 2, etc. It is easy to check that in the 
asymptotic limit v\ — > 0, Eq. (13) gives 



(n\d x \0) 

{£ n - £ y 



= -v\ 



(n\d x H\0) 
(£ n - £o) r+1 



(16) 



where all matrix elements and energies are evaluated at r = 
0. From this perturbative result we can in principle evaluate 
the leading non-adiabatic response of various observables and 
define the corresponding susceptibilities. For the purposes of 
comparing with the QAQMC approach, Eq. (16) suffices. 



B. Quasi-adiabatic quantum Monte Carlo 

The quasi-adiabatic QMC method may appear very differ- 
ent from NEQMC but has a similar underlying idea. Instead 
of imaginary time propagation with Eq. (4), we apply a simple 
operator product to evolve the initial state. We first examine 
the state propagated with the first t operators in the sequence 
P t ,i in Eq. (5), 



l^)=]l[-«(Ai)]|^) 



(17) 



and after that we will consider both symmetric expectation 
values, of the form (4> M\A\ij) m) , as well as asymmetric ex- 
pectation values Eq. (7). We are assuming that the spectrum 
of — TL is strictly positive, which is accomplished with a suit- 
able constant offset to H if needed. 



1. Linear protocols 

The coupling A can depend on the index t in an arbitrary 
way. It is convenient to define 



(18) 



where T is the overall time scale, which can be set to unity. 
The leading non-adiabatic corrections will be determined by 
the system properties and by the behavior of A(ri) near the 
point of measurement t. The most generic is the linear de- 
pendence A(t,) ps A(t) + v\(t — Ti), where v\ is related to 
the quench velocity (see below). In the end of this section 
we will briefly consider also more general nonlinear quench 
protocols. 

Our strategy to analyze Eq. (17) in the adiabatic limit will 
be the same as in the preceding subsection. We first go to the 
instantaneous basis and rewrite 

Mn)) = \i>i) = E a «( r ')l^(A i )} = 5>nl0n>- (19) 



In the instantaneous basis the discrete Schrodinger-like equa- 
tion = — - H(Ti+i)|'i/; 4 ) reads 



(20) 
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and it is instructive to compare this with Eq. (10). It is conve- 
nient to first make a transformation 



n 

j=i+ 



i 



1 {-Sl 



(21) 



This transformation does not affect the transition amplitude 
at the time of measurement t: a n = a* . Then the equation 
above becomes 



a 



i+l 



E 



n 



£ j 



Let us introduce a discrete derivative 



m l ). 



and write the Schrodinger-like equation as 

si 



n 

j=i+i 



n l \A\m l ). 



(22) 



(23) 



(24) 



In the adiabatic limit the solution of this equation is a' l n = 
S„o, i.e., the instantaneous ground state. To leading order of 
deviations from adiabaticity we find 



are very similar if £ n /£o ~ const. This can in principle al- 
ways be ensured by having a sufficiently large energy offset, 
but even with a small (or vanishing) offset we expect the ra- 
tio to be essentially constant for the range of n contributing 
significantly when the spectrum becomes gapless close to a 
quantum-critical point. If the condition indeed is properly sat- 
isfied, then from Eqs. (16) and (28) we identify the quench 
velocity as 



v x = £ A} 



(29) 



This is the main result of this section. We will confirm 
its validity explicitly in numerical studies with the QAQMC 
method in Sec. III. Since £q oc N, where N is the system 
size, we can also see that v\ oc NA\ oc N/M for a given 
total change in A over the M operators in the product. 
Let us point out that Eq. (28) can be also rewritten as 



-£ A; 



(n\<Tx\0) 
S n — So 



A A (n|&|0). 



(30) 



The first contribution here exactly matches that of Eq. (16) 
while the second term is an additional contribution corre- 
sponding to a sudden quench. 



2. Nonlinear protocols 



E 

k=0 



nf 

j=k+l c o 



(n fc |A|0 fe ), 



(25) 



where C n can be determined from the initial condition. In 
the limit of sufficiently large t the initial state is not important 
so we should have a n ~ l — > for i ^> 1, so that C n = 0. 
Therefore we find that the amplitude of the transition to the 
excited state is approximately 



«n~E 



k=0 



nf 

j=k+l c o 



(n fc |A|0 fc ). 



(26) 



Changing the summation index k to p = t — k we have 



t 

E 

P =i 



n f 



(n^ p |^|0*- p ) 



(27) 



It is clear that for large t only p<t terms contribute to the 
sum. In the extreme adiabatic limit one can thus move the 
matrix element outside of the summation and use the spectrum 
of the final Hamiltonian. In this case we find 



(n|A|0> 



So 1 — S n / Sq 

(n\K\0) 



-S n A 



S„. — Si 



(I 



e^'^M, (28) 



where A \ = \(t)—\(t—l). By comparing Eqs. (16) and (28) 
we see that near the adiabatic limit QAQMC and NEQMC 



We can extend the above result Eq. (30) to arbitrary quench 
protocols. In particular, consider 



v\ p 



r-l 



(SoY (r-l)!' 



(31) 



where r > (not necessarily an integer). For r = 1 we 
recover the linear protocol analyzed above. Then we can still 
rely on Eq. (27) but need to take into account that 



(n*-P|A"|0*- p ) w AAt_p(n*|^|0*) 



v x p 



(SoY (r - 1)! 



fy*). (32) 



Thus, we find that 



rLi 



x-rtei/fbX^I&lO*), (33) 



" (SoY(r-iy. 
where Lii_ r (x) is the Polylog function. In particular, 

1 



Li (a;) = 
Li_i(a;) 

Li_ 2 (x) 



1-x' 



(1 - a;) 2 ' 
x(x + 1) 

(1 - X) 3 ' 



(34) 
(35) 

(36) 



Under the conditions discussed above (large offset or small 
energy gap) we again have x = £ n /£o 1 and then we 
recover the continuum result using the fact that 



Lii_ r (l - e) 



(r-l)! 
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Then, indeed, 



2. Asymmetric expectation value, t 7^ M 



a' =a t 



''A 



(r - 1)1 



(-£ y(r-iy.(i-£ n /s a y 

(n*|&|0*) 



' {£n - £o) r ' 

which exactly matches Eq. (16). 



(37) 



C. Expectation values 

While asymptotically Eq. (7) gives the ground state of the 
observable A in the adiabatic limit for all values of t, the ap- 
proach to this limit as t — > 00 is qualitatively different de- 
pending on whether t is equal to M or not. More precisely 
if t = rjM where 7/ £ (0, 2) as M — » 00 we encounter two 
different asymptotic regimes for 77 ^ 1 and 77 = 1. 



7. Symmetric expectation values; t = A/ 

In this limit the expectation value of the observable A in the 
leading order of the adiabatic perturbation theory reduces to 



(A) t=M « Mvx)\A\1>(v x )), 



(38) 



where v\ rj £oAa is the imaginary time velocity identi- 
fied earlier. For generic observables not commuting with the 
Hamiltonian we find 



where 



Xaa 



= ^(0|A|n) 

n#0 



+ C.C. 



(39) 



(40) 



is the susceptibility. All energies and matrix elements are eval- 
uated at "time" t = M. 

For diagonal observables A, like the energy or energy fluc- 
tuations, we have 



(A) 



t=M 



n#0 



(£„ - £ y 



(41) 



In particular, the correction to the energy is always positive 
as it should be for any choice of wave function deviating 
from the ground state. Let us emphasize that for diagonal ob- 
servables the leading non-adiabatic response at the symmetric 
point in imaginary time coincides with that in real time, and, 
thus QAQMC or NAQMC can be used to analyze real time 
deviations from adiabaticity, as was pointed out in the case of 
NAQMC in Ref. [4]. 



It turns out that the asymptotic approach to the adiabatic 
limit is quite different for non-symmetric points t = rjM with 
^ 1, Without loss of generality we can focus on < 77 < 1. 
Then the expectation value of A is evaluated with respect to 
different eigenstates; 



(A)t 



(i/> L \A\t/> R ) 
(iPl\iPr) 



(42) 



where 



\^r) = Hn(\ k )\i> ), 

fc=l 
t+1 

Ml) = [ n(x k )\^ } 



(43) 



k=M 



Note that the overlap (ipL\ipn) is independent of t by con- 
struction. It is real for symmetric protocols \ m -k = which 
we consider. 

It is easy to see that for diagonal observables we obtain 
leading asymptotic as in Eq. (41) but with the opposite sign 
in the second term 



\(n\dx\0)\ 2 

(£n ~ £o)< 



(n\A\n) 



(44) 



In particular, the leading correction to the ground state energy 
is negative when t deviates sufficiently from the symmetric 
point. There is no contradiction here since the left and right 
states are different (i.e., we are not evaluating a true expecta- 
tion value and there is no variational principle). Since the cor- 
rection up to the sign exactly matches the real time result we 
can still use the non-symmetric expectation value for diagonal 
observables to extract the real time non-adiabatic response. 
For t —> M the sign of the correction should change, to con- 
nect smoothly to the variational t = M expectation value. 

As in the symmetric case, using the APT discussed in the 
previous section the results derived here easily extend to other 
values of the exponent r. 



3. The metric tensor and Berry curvature 

If A = —d^H, then the susceptibility Eq. (40) reduces to 
the p;A component of the metric tensor, 4,7 which, thus, can be 
readily extracted using the QAQMC algorithm. In particular, 
the diagonal components of the metric tensor define the more 
familiar fidelity susceptibility. 

Next, we observe that for t sufficiently different from M, 
the approach to the ground state in the left function ipL in 
Eq. (43) effectively corresponds to a change in sign of the 
velocity, and, thus, we find 



(A}t 



(if,(-v)\A\j>(v)) 

W-t>)|iK«)) 



(45) 
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We can use the results of the previous section to find that for 
off-diagonal observables 

(A) t « (A) - ivxXax, (46) 

Xaa = * >,<0|^4.|n>— — - c.c. (47) 

Based on this result we conclude that the leading non- 
adiabatic correction is imaginary and coincides, up to the 
factor of imaginary i, with the real-time non-adiabatic 
correction. 7 In particular, for A = —d^H the susceptibility 
x" AX = XuA ls proportional to the Berry curvature. We will 
discuss this interesting fact further in Ref. [7]. 



III. RESULTS 

As a demonstration of the utility of QAQMC and the behav- 
iors derived in the previous section we here study the TFIM, 
defined by the Hamiltonian 

U = -sY J <7?<7?-(l- s )5>f, (48) 

(i.j) 

where are nearest-neighbor sites, and a z and a x are 

Pauli matrices. Here s plays the role of the tuning param- 
eter, which in the simulations reported below will vary be- 
tween (where the ground state is trivial) to a value exceed- 
ing the quantum-critical point; s c = 1/2 in a ID chain and 
s c « 0.247 in the 2D square lattice. 10 

We work in the standard basis of eigenstates of all Sf. The 
simulation algorithm samples strings of 2M diagonal and off- 
diagonal terms in Eq. (48), in a way very similar to the T > 
stochastic series expansion (SSE) method, which has been dis- 
cussed in detail in the case of the TFIM in Ref. [11]. The mod- 
ifications for the QAQMC primarily concern the sampling of 
the initial state, here ^(O)) = ]X | fj + |;), which essen- 
tially amounts to a particular boundary condition replacing 
the periodic boundaries in finite-temperature simulations. An 
SSE-like scheme with such modified boundaries were also 
implemented for the NEQMC method in Ref. [4], and re- 
cently also in a study of combinatorial optimization problems 
in Ref. [12]. We here follow the same scheme, using cluster 
updates in which clusters can be terminated at the boundaries. 
The implementation for the product with varying coupling s 
is even simpler than SSE or NEQMC, with the fixed-length 
product replacing the series expansion of Eq. (4). The changes 
relative to Refs. [4] and [11] are straightforward and we there- 
fore do not discuss the sampling scheme further here. 



A. Quantum-critical dynamic scaling 

The idea of dynamic scaling at a critical point dates back to 
Kibble and Zurek for quenches (or ramps, since the parameter 
does not change suddenly, but linearly as a function of time) 



through classical phase transitions. • These ideas were later 
extended also to quantum systems. 15-18 The basic notion is 
that the system has a relaxation time t Te u an d if some parame- 
ter (here a parameter of the Hamiltonian) is changed such that 
a critical point is approached, the system can stay adiabatic (or 
in equilibrium) only if the remaining time t to reach the criti- 
cal point is much larger than the relaxation time; t S> t Te \. In 
general one expects t Ie \ ~ £ 2 ~ e" 2I/ where £ is the correla- 
tion length, v the exponent governing its divergence with the 
distance e to the critical point, and z the dynamic exponent. 
For a system of finite size (length) L, £ is maximally of order 
L and, thus, for a linear quench the critical velocity v CT ft sepa- 
rating slow and fast power-law quenches according to Eq. (15) 
should heuristically be given by v cr n ~ L~( z+1 ' v \ and for a 
power-law quench with exponent r according to Eq. (15) this 
generalizes to 9 

v crit ~L-(* r +i\ (49) 

One then also expects a generalized finite-size scaling form 
for singular quantities A; 

A(L, e) = L K f(eL xl \ vL zr+1/ »), (50) 

where k characterizes the leading size-dependence at the crit- 
ical point of the quantity considered. For v — > Eq. (50) 
reduces to the standard equilibrium finite-size scaling hypoth- 
esis. This scaling was recently tested in different systems, 
both quantum 9 ' 19,20 and classical 21 . 

The above expression Eq. (50) combined with the product- 
evolution Eq. (5) allows us to study a phase transition based 
on different combinations of scaling in the system size and 
the velocity in non-equilibrium setups. For example, if one 
wants to find the critical point for the phase transition and the 
exponent v is known, one can carry out the evolution under 
the critical-velocity condition: 

vL z+ »=c, (51) 

where c is a constant. In this paper we focus on linear quench 
protocols and set r = 1 henceforth. As we discussed in 
Sec. II B, the QAQMC method applied to a system of size 
(volume) N based on evolution with M operators in the se- 
quence and change A a between each successive operator cor- 
responds to a velocity v oc NA\ oc N/M, with the prefactor 
depending on the ground state energy (at the critical point). 
The exact prefactor will not be important for the calculations 
reported below, and for convenience in this section we define 



where Sf is the final value of the parameter s in Eq. (48) over 
the evolution (which is also the total change in s, since we 
start with the eigenstate at s = 0). The critical product-length 
M is, thus, given by 

M = -NL z+1/lJ = -L d+Z+ V v , (53) 
c c 

where we have also for simplicity absorbed Sf into c. 
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Using an arbitrary c of order 1 in Eq. (51), the critical point 
s c can be obtained based on a scaling function with the single 
argument eL 1 ^ in Eq. (50). We will test this approach here, in 
Sees. IIIC and HID. and later, in Sec. HIE, we will show that 
exact knowledge of the exponents in Eq. (51) are actually not 
needed. First we discuss the quantities we consider in these 
studies. 



B. Quantities studied 

We will focus our studies here on the squared z-component 
magnetization (order parameter), 

i / N \ 2 

<=(M^v >' (54) 

We can also define a susceptibility-like quantity (which we 
will henceforth refer to as the susceptibility) measuring the 
magnetization fluctuations: 

X = N((m 2 z ) - (\m z \) 2 ). (55) 

Here both terms have the same critical size-scaling as the 
equal-time correlation function; 

<m z ) 2 ~ (\m z \f ~ L^ d +^+v) i (56) 

where d is the spatial dimensionality. The prefactors for the 
two quantities are different, however, a divergent peak re- 
mains in Eq. (55) at the transition. Away from the critical 
point X — >■ with increasing system size. 

To clarify our use of x, we point out that we could also 
just study the scaling of (to 2 ,), but the peak produced when 
subtracting off the second term in Eq. (55) is helpful in the 
scaling analysis. According to Eq. (50) and using z = 1 in 
Eq. (56), the full scaling behavior of the fluctuation around 
the critical point should follow the form 

X~£ 1 - , '/((.s- Sc )£ 1 / l >£ 1+ ^), (57) 

for any dimensionality d. 

We should point out here that the true thermodynamic sus- 
ceptibility based on the Kubo formula 22 (imaginary-time in- 
tegral) yields a stronger divergence L d ~ L 2 ~ n . This quan- 
tity is, however, more difficult to study with the QAQMC al- 
gorithm, because, unlike in standard finite-T QMC methods, 
the time integration cannot simply be carried out within the 
space of time-evolving Hamiltonians in Eq. (5) and Eq. (7). 
The standard Feynman-Suzuki correspondence between the 
d-dimensional quantum and -dimensional classical sys- 

tems is not realized in our scheme. The configuration space of 
time-evolving Hamiltonians builds in the relaxation time, t re \, 
in a different way, not just in terms of equilibrium fluctuations 
in the time direction, but in terms of evolution as a function of 
a time-dependent parameter. 

A useful quantity to consider for extracting the critical point 
is the Binder cumulant; 23 , 

3 / 1 (mi) \ 
U=- (1--^— *4 J. (58) 

2 V 3 (TO 2 ) 2 / 



For a continuous phase transition, U converges to a step func- 
tion as L —> oo. The standard way to analyze this quantity 
for finite L is to graph it versus the scaling argument eL x / v 
for different L and extract crossing points, which approach 
the critical point with increasing L. Normally this is done in 
the equilibrium, either by taking the limit of the temperature 
T — > for each L first, or by fixing (3 = 1/T cx L z if z is 
known. Here the latter condition is replaced by Eq. (51), but, 
as we will discuss further below, the condition can be relaxed 
and the exponents do not have to be known accurately a priori. 
Our approach can also be used to determine the exponents, ei- 
ther in a combined procedure of simultaneously determining 
the critical point and the exponents, or with a simpler analysis 
after first determining the critical point. 

We have up until now only considered calculations of 
equal-time observables, but in principle it is also possible 
and interesting to study correlations in the evolution direction, 
which also can be used to define susceptibilities. 

In the following we will illustrate various scaling proce- 
dures using results for the ID and 2D TFIMs. The dynamic 
exponent z = 1 is known for for both cases, and in the ID 
case all the exponents are rigorously known since they coin- 
cide with those of the classical 2D Ising model. For the 2D 
TFIM the exponents are know rather accurately based on nu- 
merics for the 3D classical model. 



C. ID transverse-field Ising model. 

The ID TFIM provides a rigorous testing ground for the 
new algorithm and scaling procedures since it can be solved 
exactly. 24 The critical point corresponds to the ratio between 
the transverse field and the spin-spin coupling equaling 1, i.e., 
s = 1/2 in the Hamiltonian Eq. (48). The critical exponents, 
known through the mapping to the 2D Ising model, 25 are v = 
1 and rj = 1/4. 

The results presented here were obtained in simulations 
with the parameter s spanning the range [0, S{] with S{ = 0.6, 
i.e., going from the trivial ground state of the field term to 
well above the critical point. Fig. 1 shows examples of results 
for the susceptibility and the Binder cumulant. The operator- 
sequence length M, Eq. (5), was scaled with the system size 
in order to stay at the critical velocity according to Eq. (53). 
We emphasize again that a single run produces a full curve 
within the s-range used. In order to focus on the behavior 
close to criticality we have left out the results for small s in 
Fig. 1. Since M is very large (up to w 10 6 for the largest 
L in the cases shown in the figure), we also do not compute 
expectation values for each t in Eq. (7), but typically spacing 
measurements by oc TV operators. 

Extracting Binder curve-crossings using system-size pairs 
L and L + 4, with L = 4, 8, 12, . . . 60, and extrapolating the 
results to L —> oo, we find s c — 0.49984(16), as illustrated 
in Fig. (2). Thus, the procedure produces results in full agree- 
ment with the known critical point. 

The dynamical scaling of the susceptibility is illustrated in 
Fig. 3. Here there are no adjustable parameters at all, since 
all exponents and the critical coupling are known (and we use 
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FIG. 1: (Color online) Results of typical QAQMC runs for the ID 
TFIM, Eq. (48). The binder cumulant Eq. (58) (bottom panel) and the 
susceptibility \ Eq. (55) (top panel) are graphed versus s for several 
system sizes L. In these simulations, which spanned the range s £ 
[0, 0.6], the length of the index sequence was of the form Eq. (53), 
i.e., with the exponents applicable in this case M — L 3 /c with the 
arbitrary constant chosen to be c = 4 3 /240. 



the exact critical coupling s c = 1/2, although the numerical 
result extracted below is very close to this value and produces 
an almost identical scaling collapse). While some deviations 
from a common scaling function are seen for the smaller sys- 
tems and far away from the scaled critical point (s — s c )L, the 
results for larger sizes and close to the peak rapidly approach a 
common scaling function. This behavior confirms in practice 
our discussion of the definition of the velocity and the ability 
of the QAQMC method to correctly take into account at least 
the first corrections to the adiabatic evolution. 



D. 2D transverse-field Ising model 

The 2D transverse-field Ising model provides a more seri- 
ous test for our algorithm since it cannot be solved exactly. 
Among many previous numerical studies, 10,26 ' 27 Ref. [10] ar- 
guably has the highest precision so far for the value of the 
critical coupling ratio. Exact diagonalization was there car- 
ried out for up to 6 x 6 lattice size. In terms of the critical field 
h c = 1 — s in units of the coupling J = s, the critical point 
was determined to h c /J = 1/0.32841(2) = 3.04497(18), 
where the error bar reflects estimated uncertainties in finite- 
size extrapolations. Results based on QMC simulations 26,27 
are in agreement with this value, but the statistical errors are 
larger than the above extrapolation uncertainty. One might 




0.04 0.05 
1/L 

FIG. 2: (Color online) Results of a Binder-crossing scaling analysis 
of the ID TFIM data in Fig. 1 (including also other system sizes not 
shown there). Crossing points were extracted based on system sizes 
L and L + 4, with L = 4,8, ... , 60. The curve is a fit to the form 
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0.49984(16) andb = 1.6(1). 



worry that the system sizes L < 6 are very small and the 
extrapolations may not reflect the true asymptotic L — > oo 
size behavior. However, the data points do follow functional 
forms expected based on the corresponding low-energy field 
theory, and there is therefore no a-priory reason to question 
the results. It is still useful to try to reach similar or higher 
precision with other approaches, as we will do here with the 
QAQMC method combined with dynamic scaling. 

In this case we simulate the linear quench in the window 
of s G [0, 0.3], which contains the previous estimates for the 
critical value s c « 0.247 as discussed above. Although we 
could also carry out an independent scaling analysis to extract 
the critical exponents, we here choose to just use their val- 
ues based on previous work on the classical 3D Ising model; 
\jv « 1.59, and rj « 0.036. 28 Our goal here is to extract 
a high-precision estimate of the critical coupling, and at the 
same time to further test the ability of QAQMC to capture the 
correct critical scaling behavior. We again scale M with L 
according to Eq. (53), with the constant c = 4 4 59 /32. 

As in the ID case, we extract Binder-cumulant crossing 
points based on linear system sizes L and L + 4 with L = 
4, 8, ... , 56. Fig. 4 shows the results versus 1/L along with 
a fit to a power-law correction for s c (L). Extrapolating to 
infinite size gives s c = 0.247241(7), which corresponds to a 
critical field (in unit of J) h c / J = 3.04463(12). This is in rea- 
sonable good agreement with the value obtained in Ref. [10] 
and quoted above, with our (statistical) error bar being some- 
what smaller. To our knowledge, this is the most precise value 
for the critical coupling of this model obtained to date. We 
emphasize that we here relied on the non-equilibrium scaling 
ansatz to extract the equilibrium critical point. Allowing for 
deviations from adiabaticity in a controlled way and utilizing 
the advantage of the QAQMC algorithm allowed us to extract 
observables in the whole range of couplings in a single run. 
This requires considerably less computational resources than 
standard equilibrium simulations, which must be repeated for 
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FIG. 3: (Color online) Scaled susceptibility of the ID TFIM. The 
axes have been scaled according to the form (57) with the second 
argument constant and using the exact critical point s c — 1/2. The 
results are shown on two different scales to make visible deviations 
(due to subleading size and velocity corrections) from the common 
scaling function far away from criticality as well as the good data 
collapse close to the critical point. 
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FIG. 4: (Color online) Binder crossings for the 2D TFIM extracted 
using L and L + 4 systems with L = 4, 8, ... , 56. The crossing 
points have been fitted to the standard form s c (L) = s c + a/L b , for 
which the optimal values are s c = 0.247241(7) and b = 4.2(4). 
The results are shown on two different scales to illustrate large devi- 
ations from the fitted form for the smaller systems, followed by rapid 
convergence for larger sizes. 



several different couplings in order to carry out the crossing- 
point analysis. 

Fig. 5 shows the susceptibility scaled according to the be- 
havior expected with Eq. (50) when the second argument is 
held constant. As in the ID case, the data converge rapidly 
with increasing size toward a common scaling function in the 
neighborhood of the transition point, again confirming the cor- 
rect quasi-adiabatic nature of the QAQMC method. 



E. Further tests 

The results discussed in the preceding subsections were ob- 
tained with the KZ velocity condition Eq. (51), applied in the 
form of Eq. (53) tailored to the QAQMC approach, with spe- 
cific values for the constant c. In principle the constant is 
arbitrary, but the non-universal details of the scaling behavior 
depend on it. This is in analogy with a dependence on the 
shape, e.g., an aspect ratio, of a system in equilibrium simu- 
lations at finite temperature, or to the way the inverse temper- 
ature (3 = 1 /T is scaled as aL z with arbitrary a in studies of 
quantum phase transitions (as an alternative to taking the limit 
j3 — > oo for each lattice size). The critical point and the criti- 
cal exponents should not depend on the choices of such shape 
factors or limiting procedures. 

To extract the critical coupling, in the preceding subsec- 
tions we fixed the exponents v and z at their (approximately) 
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FIG. 5: (Color online) Scaled susceptibility of the 2D TFIM, based 
on Eq. (57) with a constant second argument. Here we have used 
1/v — 1.59 and r] = 0.036 for the classical 3D Ising model 28 . 
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known values, and one may at first sight assume that it is nec- 
essary to use their correct values. It is certainly some times 
convenient to do so, in order to set the second argument of the 
scaling function Eq. (50) to a constant and, thus, obtain a sim- 
pler scaling function depending on a single argument. How- 
ever, one can study critical properties based on the scaling 
approach discussed above as long as the velocity approaches 
zero as the system size increases. This observation can be im- 
portant in cases where the critical exponents are not known 
and one would like to obtain an accurate estimate of the crit- 
ical coupling before carrying out a scaling analysis to study 
exponents. We will test this in practice here. As we will dis- 
cuss further below, one should use a different power k in the 
scaling ansatz Eq. (50) if the velocity is brought to zero slower 
than the critical form. 

In cases where we use the "wrong" values of the exponents, 
we formally replace z + 1/v by a free parameter a, 



L- a /c, 



(59) 



and the corresponding substitution in Eq. (53). To understand 
the scaling of the observables for arbitrary a we return to 
the general scaling form given by Eq. (50). In the case of 
the Binder cumulant and for linear quench protocol, this form 
reads 



U = f(ts-s c )L l / v ,vL*+ 1 l»). 



(60) 



As we pointed out above, when the velocity scales exactly 
as L _ ( 2+1 / ,y ) the dependence on the second argument in the 
scaling function drops out and we can find the crossing point 
in a standard way as we did in Figs. 2 and 4. Suppose that we 
do not know the exponents v and z a priory and instead scale 
v as in Eq. (59). Then there are three possible situations: (i) 
a = z + 1/v, (ii) a > z+l/v and (iii) a < z + 1/v, where 
we already have analyzed scenario (i). In scenario (ii), where 
velocity scales to zero faster than the critical KZ velocity, 
the second argument of the scaling function L z+1 / V / L a ap- 
proaches zero as the system size increases and, thus, the scal- 
ing function effectively approaches the equilibrium velocity- 
independent form. We can then extract the crossing point as in 
the first scenario, and this gives the correct critical coupling in 
the limit of large system sizes. Finally, in case (iii) the velocity 
scales zero slower than the critical KZ value and the second 
argument in Eq. (60) diverges, which implies that the system 
enters a strongly non-equilibrium regime. This scenario ef- 
fectively corresponds to taking the thermodynamic limit first 
and the adiabatic limit second. Then, if the system is initially 
on the disordered side of the transition, the Binder cumulant 
vanishes in the thermodynamic limit. At finite but large sys- 
tem sizes its approach to zero should be given by the standard 
Gaussian theory: 



(61) 



Combining this with the scaling ansatz Eq. (60) we find that 
for a < z + l/v the expected asymptotic of the Binder cumu- 
lant is 
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FIG. 6: (Color online) Critical-point estimates based on curve cross- 
ings of appropriately scaled quantities for scenarios (i), (ii), and (iii) 
discussed in the text. The Binder cumulant (lower panel) and the 
squared magnetization (upper panel) give estimates s c (L) and s' c (L), 
respectively, based on system sizes L and L + 4. The red and blue 
curves correspond to runs in which the velocity was kept at the crit- 
ical value, scenario (i), but with different constants of proportional- 
ity c in Eq. (53); ci = 4 4 ' 59 /32 and c 2 = 4 4 ' 59 /48. The yellow 
curves were obtained with the velocity decreasing faster than v c with 
L, scenario (ii), with the proportionality constant c;; = 4 /32. The 
green and pink curves correspond to cases where the velocity is sub- 
critical, scenario (iii), with constants C4 = 4 4,2 /32, C5 = 4 4 /32. 
In all cases, power-law corrections were fit in order to extrapolate to 
infinite size (with small sizes excluded until statistically sound fits 
were obtained). 



where / is some other velocity independent scaling function. 
Thus we can find the correct transition point by finding cross- 
ing points of U L d v d ^ z+1 / y \ Similar considerations apply to 
the ordered side of the transition, where the Binder cumulant 
approaches one as the inverse volume. 

The three cases are illustrated in the lower panel of Fig. 6, 
which shows Binder-cumulant crossings extracted from ap- 
propriately scaled data in cases (i), (ii), and (iii) above. Ad- 
ditionally, to illustrate the insensitivity to the choice of the 
constant c in the scaled sequence length in Eq. (53), results 
based on two different constants are shown for case (i). In 
all cases, the extrapolated critical couplings agree with each 
other to within statistical errors. Note that, one the one hand, 
if the exponent a gets very large, then the time of simulations, 
which scales as M, rapidly increases with the system size and 
the algorithm becomes inefficient. On the other hand, if a 
is very small our results indicate that the size-dependence is 
larger and it is more difficult to carry out the extrapolation to 
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infinite size. The optimal value of a should be as close as pos- 
sible to the critical KZ power, but to be on the safe side when 
scaling according to the standard KZ critical form, case (i), 
one may choose a somewhat larger value, since the subcritical 
velocity case (ii) has the same scaling form. 

Next we illustrate how the same idea works in the case of 
the order parameter. Around the critical point (s c , v c ritX the 
squared magnetization [see Eq. (54)] can be written as 

ml = L-^l»f{(s- Sc )l}l\vL* +1 l v ). (63) 

As in the previous discussion we scale v ~ L~ a and depend- 
ing on the exponent a there are two different asymptotics of 
the scaling function. For a > z + 1/v the second argument 
vanishes or approaches constant so we effectively get the equi- 
librium scaling 

m 2 = L -W/u f ( {s _ Sc)L X/v} (M) 

If, conversely, a < z + 1 jv then on the disordered side of the 
transition to 2 , scales as L~ d . This immediately determines the 
asymptotic of the scaling function in Eq. (63): 

m 2 = L - dv S2§{$=£ _ Sc)L i/vy (65) 

Equation (65) can be used in the same way as the Binder 
cumulant to extrapolate the critical point, using the standard 
form s' c (L) = s' c + a/L b for the rescaled to 2 . As shown in the 
top panel of Fig. (6), after rescaling the order parameter and 
extrapolating the crossing points between the appropriately 
rescaled to 2 curves to the thermodynamics limit, all curves, 
obtained from below or above the adiabatic limit Eq. (49), 
converge to the same value s' c « 0.247. This approach also 
suggests a way to determine the transition point in experiment, 
since one can sweep through the critical point at different ve- 
locities, the crossing point can then be extracted through the 
measurement of the order parameter. It's also worth mention- 
ing that since one can extrapolate the critical point indepen- 
dently without knowing the actual exponent v prior to the 
simulation, an optimization procedure can be carried out to 
determine the exponents posterior to the simulation. 29 

For completeness we also briefly discuss the role of the fi- 
nal point Sf of the evolution. Fig. 7 shows ID results for the 
squared magnetization and susceptibility obtained for a range 
of final points above the critical value. Here the velocity was 
kept constant for all the cases. The values of the computed 
quantities at some fixed s, e.g., at s c , show a weak dependence 
on Sf for the lowest-Sf runs. The deviations are caused by con- 
tributions of order v 2 and higher, which are non-universal as 
discussed in Sec. II B. For very high velocities the dependence 
on Sf can much more be dramatic than in Fig. 7, but this is not 
the regime in which the QAQMC should be applied to study 
universal physics. 

IV. SUMMARY AND DISCUSSION 

We have presented a nonequilibrium QAQMC approach to 
study quantum dynamics, with a simple product of operators 
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FIG. 7: (Color online) Squared magnetization (lower panel) and sus- 
ceptibility (upper panel) versus s of the ID TFIM with L — 24. In 
these runs, different curves correspond to different end points Sf of 
the evolution, with the velocity v oc SfN/M kept constant. 



with evolving coupling replacing the standard Schrodinger 
time evolution. We showed that this approach captures the 
leading non-adiabatic corrections to the adiabatic limit, both 
by analytical calculations based on the APT and by explicit 
simulations of quantum-critical systems with the QAQMC al- 
gorithm. The simulation results obey expected generalized 
dynamic scaling with known static and dynamic critical expo- 
nents. We also extended the scaling formalism beyond results 
obtained previously in Ref. [4]. We analyzed the leading non- 
adiabatic corrections within this method and showed that they 
can be used to extract various non-equal time correlation func- 
tions, in particular the Berry curvature and the components of 
the metric tensor. A clear advantage of the QAQMC approach 
is that one can access the whole range of couplings in a sin- 
gle run. Being a simple modification of projector QMC, the 
QAQMC method is applicable to the same class of models as 
this conventional class of QMC schemes — essentially models 
for which "sign problems" can be avoided. 

As an illustration of the utility of QAQMC, we applied the 
algorithm and the scaling procedures to the ID and 2D TFIMs. 
The expected scaling behaviors are observed very clearly. In 
the ID case we extracted a critical coupling in full agreement 
with the known value, and in 2D we obtained an estimate 
with unprecedented (to our knowledge) precision (small error 
bars); (h/J) c = 3.04463(12). 

The QAQMC approach bears some similarities to pre- 
vious implementations of quantum annealing within QMC 
algorithms. 30,31 However, the previous works have mainly 
considered standard equilibrium QMC approaches in which 
some system parameter is changed as a function of the sim- 
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illation time. This evolution is not directly related to true 
quantum dynamics (and, thus, is not really quantum anneal- 
ing), but is dependent on the particular method used to up- 
date the configurations. In contrast, in our scheme, as in the 
NEQMC method introduced in Ref. [4], the evolution takes 
place within the individual configurations, and there is a di- 
rect relationship to true Schrodinger evolution in imaginary 
time. In Green's function (GF) QMC simulations the grad- 
ual change of a system parameter with the simulation time 
is rather closely related to the QAQMC scheme (since also 
there one applies a series of terms of the Hamiltonian to a 
state), with the difference being QAQMC uses true impor- 
tance sampling of configurations, with no need for guiding 
wave functions and no potential problems related to mixed es- 
timators. Our asymmetric expectation values could be consid- 
ered as a kind of mixed estimator as well, but we have com- 
pletely characterized them within the APT. In addition, the 
previous uses of GFQMC with time-evolving Hamiltonians 
have, to our knowledge, never addressed the exact meaning 
of the velocity of the parameter evolution. The correct defini- 
tion of the velocity is of paramount importance when applying 
quantum-critical scaling methods, as we have discussed here. 
We have here computed the velocity within within APT for the 
QAQMC scheme. The same relationship with Schrodeinger 
dynamics may possibly hold for GFQMC as well, but, we 
have not applied the APT to this case and it is therefore not 
yet clear whether GFQMC can capture correctly the same 
universal non-equilibrium susceptibilities as the QAQMC and 
NEQMC methods. We expect QAQMC to be superior to time- 
evolving GFQMC, because of its better control over measured 
symmetric and asymmetric expectation values and fully real- 
ized importance sampling. 

We also stress that we have here not focused on optimiza- 
tion. Previous works on quantum annealing within QMC 
schemes have typically focused on their abilities to optimize 
difficult classical problems. While the QAQMC may poten- 
tially also offer some opportunities in this direction, our pri- 
mary interest in the method is to use it to extract challenging 
dynamical information under various circumstances. 

The QAQMC and NEQMC methods provide correct real- 
izations of quantum annealing in imaginary time. Besides 
their ability to study dynamic scaling, with exponents iden- 
tical to those in real-time Schrodeinger dynamics, 4 it will be 
interesting to explore what other aspects of real-time dynam- 
ics can be extracted with these methods. In particular, their 
applicability to quantum glasses, of interest in the context of 
quantum adiabatic computing 12 as well as in condensed mat- 
ter physics, deserves further studies. 

The ability of the QAQMC to produce results for a whole 
evolution path in a single run can in principle also be carried 
over to the conventional Schrodeinger imaginary-time evolu- 
tion with U{t) in Eq. (4). By "slicing" the time evolution into 
K successive evolutions over a time-segment A T ; 
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where r n = nA r , one can evaluate matrix elements anal- 
ogous to Eq. (7) by inserting the operator of interest at 



FIG. 8: (Color online) One-way evolution s £ [0, 1] with QAQMC 
for the ID TFIM. Lower panel: The susceptibility Eq. (55). Upper 
panel: the rescaled susceptibility Eq. (57). Each full curve corre- 
sponding to a given chain length L was obtained in a single run. The 
velocity Eq. (52) was held constant at 4 3 /80. 



any point within the product of time-slice operators in 
(*(A )|C/*(r)C/(r)|*(A )). In this case, the symmetric ex- 
pectation value, evaluated at the mid-point, is identical to 
the NEQMC method, 4 and the asymmetric expectation values 
will exhibit properties similar to those discussed in Sec. II C 2. 
We have not yet explored this approach, and it is not clear 
whether it would have any other advantage besides the exact 
reduction to Schrodinger dynamics of the symmetric expecta- 
tion values. In practice the simulations will be more complex 
than the QAQMC approach because of the need to sample in- 
tegrals, but not much more so than the NEQMC method. 

Finally, we point out that in principle one can also carry 
out a one-way evolutions with the QAQMC algorithm. In- 
stead of starting with the A = Ao eigenstate at both (■)/>£ | and 
\ipn) and then projecting them to the A = Xm eigenstate us- 
ing two sequences of the form Eq. (5), one can make (iPl\ 
correspond to Ao and let it evolve to \ipR) corresponding to 
Am with only a single operator sequence of length M. In 
the case of the TFIM Eq. (48), the obvious choice is then to 
evolve from s = to s = 1 (the classical Ising model), so 
that both edge states are trivial. All our conclusions regarding 
the definition of the velocity and applicability of scaling form 
remain valid in this one-way QAQMC. Results demonstrat- 
ing this in the case of the ID TFIM are sown in Fig. 8. We 
anticipate that this approach may be better than the two-way 
evolution in some cases, but we have not yet compared the 
two approaches extensively. 
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